"""
Material Properties Framework

Implementation of temperature-dependent and pressure-dependent material properties
for multi-physics CFD simulations. Provides comprehensive property databases for
common fluids and gases, with support for custom property correlations.

Key features:
- Temperature-dependent viscosity, thermal conductivity, and density
- Pressure-dependent properties for compressible flows
- Built-in property databases for common fluids
- Custom property correlation support
- Non-Newtonian fluid models
- Mixture property calculations
- Property validation and bounds checking

Supported Fluids:
- Air (ideal gas)
- Water (liquid)
- Steam (gas phase water)
- Common engineering fluids
- Custom fluid definitions

Property Models:
- Sutherland's law for viscosity
- Power law correlations
- Polynomial correlations
- Piecewise linear interpolation
- Ideal gas law for density
- Custom functional forms
"""

using LinearAlgebra

# Abstract types for extensibility
abstract type FluidProperties end
abstract type PropertyModel end

"""
    ConstantProperties <: FluidProperties

Constant fluid properties (temperature and pressure independent).

# Fields
- `density::Float64`: Density ρ (kg/m³)
- `viscosity::Float64`: Dynamic viscosity μ (Pa·s)
- `thermal_conductivity::Float64`: Thermal conductivity k (W/m·K)
- `specific_heat::Float64`: Specific heat capacity cₚ (J/kg·K)
- `thermal_expansion::Float64`: Coefficient of thermal expansion β (1/K)
- `compressibility::Float64`: Isothermal compressibility κ (1/Pa)
"""
Base.@kwdef struct ConstantProperties <: FluidProperties
    density::Float64 = 1.225                      # kg/m³ (air at STP)
    viscosity::Float64 = 1.81e-5                  # Pa·s (air at STP)
    thermal_conductivity::Float64 = 0.025         # W/m·K (air at STP)
    specific_heat::Float64 = 1007.0               # J/kg·K (air)
    thermal_expansion::Float64 = 3.67e-3          # 1/K (air at STP, ≈ 1/T)
    compressibility::Float64 = 1e-5               # 1/Pa (typical for liquids)
end

"""
    TemperatureDependentProperties <: FluidProperties

Temperature-dependent fluid properties with various correlation models.

# Fields
- `viscosity_model::PropertyModel`: Viscosity correlation model
- `thermal_conductivity_model::PropertyModel`: Thermal conductivity correlation
- `density_model::PropertyModel`: Density correlation model
- `specific_heat_model::PropertyModel`: Specific heat correlation
- `reference_conditions::NamedTuple`: Reference state (T₀, P₀, ρ₀, μ₀, etc.)
- `temperature_range::Tuple{Float64,Float64}`: Valid temperature range (K)
- `pressure_range::Tuple{Float64,Float64}`: Valid pressure range (Pa)
"""
Base.@kwdef struct TemperatureDependentProperties <: FluidProperties
    viscosity_model::PropertyModel
    thermal_conductivity_model::PropertyModel
    density_model::PropertyModel
    specific_heat_model::PropertyModel
    reference_conditions::NamedTuple
    temperature_range::Tuple{Float64,Float64} = (200.0, 2000.0)  # K
    pressure_range::Tuple{Float64,Float64} = (1e3, 1e7)         # Pa
end

"""
    NonNewtonianProperties <: FluidProperties

Non-Newtonian fluid properties with shear-rate dependent viscosity.

# Fields
- `base_properties::FluidProperties`: Base fluid properties
- `rheology_model::Symbol`: Rheological model (:power_law, :bingham, :herschel_bulkley)
- `rheology_parameters::Dict{Symbol,Float64}`: Model-specific parameters
- `temperature_dependence::Bool`: Whether rheological parameters depend on temperature
"""
Base.@kwdef struct NonNewtonianProperties <: FluidProperties
    base_properties::FluidProperties
    rheology_model::Symbol = :power_law
    rheology_parameters::Dict{Symbol,Float64} = Dict{Symbol,Float64}()
    temperature_dependence::Bool = false
end

# Property correlation models

"""
    SutherlandModel <: PropertyModel

Sutherland's law for temperature-dependent viscosity:
μ(T) = μ₀ * (T/T₀)^(3/2) * (T₀ + S)/(T + S)

Where S is the Sutherland constant.

# Fields
- `reference_viscosity::Float64`: Reference viscosity μ₀ (Pa·s)
- `reference_temperature::Float64`: Reference temperature T₀ (K)
- `sutherland_constant::Float64`: Sutherland constant S (K)
"""
Base.@kwdef struct SutherlandModel <: PropertyModel
    reference_viscosity::Float64       # μ₀ (Pa·s)
    reference_temperature::Float64     # T₀ (K)
    sutherland_constant::Float64       # S (K)
end

"""
    PowerLawModel <: PropertyModel

Power law correlation: Property(T) = A * T^n

# Fields
- `coefficient::Float64`: Coefficient A
- `exponent::Float64`: Exponent n
"""
Base.@kwdef struct PowerLawModel <: PropertyModel
    coefficient::Float64  # A
    exponent::Float64      # n
end

"""
    PolynomialModel <: PropertyModel

Polynomial correlation: Property(T) = Σ(aᵢ * T^i)

# Fields
- `coefficients::Vector{Float64}`: Polynomial coefficients [a₀, a₁, a₂, ...]
"""
Base.@kwdef struct PolynomialModel <: PropertyModel
    coefficients::Vector{Float64}
end

"""
    PiecewiseLinearModel <: PropertyModel

Piecewise linear interpolation from tabulated data.

# Fields
- `temperature_points::Vector{Float64}`: Temperature points (K)
- `property_values::Vector{Float64}`: Property values at temperature points
"""
Base.@kwdef struct PiecewiseLinearModel <: PropertyModel
    temperature_points::Vector{Float64}
    property_values::Vector{Float64}
end

"""
    IdealGasModel <: PropertyModel

Ideal gas law for density: ρ = P*M/(R*T)

# Fields
- `molar_mass::Float64`: Molar mass M (kg/mol)
- `gas_constant::Float64`: Universal gas constant R (8.314 J/mol·K)
"""
Base.@kwdef struct IdealGasModel <: PropertyModel
    molar_mass::Float64 = 0.02897      # kg/mol (air)
    gas_constant::Float64 = 8.314      # J/mol·K
end

"""
    CustomFunctionModel <: PropertyModel

Custom user-defined property function.

# Fields
- `property_function::Function`: Function f(T, P=1e5) that returns property value
- `parameter_dict::Dict{String,Any}`: Additional parameters for the function
"""
Base.@kwdef struct CustomFunctionModel <: PropertyModel
    property_function::Function
    parameter_dict::Dict{String,Any} = Dict{String,Any}()
end

# Pre-defined fluid property databases

"""
    create_air_properties(; model=:temperature_dependent) -> FluidProperties

Create air property database with temperature-dependent correlations.

# Arguments
- `model::Symbol`: Property model type (:constant, :temperature_dependent, :ideal_gas)

# Returns
- `FluidProperties`: Air properties suitable for atmospheric conditions
"""
function create_air_properties(; model::Symbol=:temperature_dependent)
    if model == :constant
        return ConstantProperties(
            density = 1.225,           # kg/m³ at STP
            viscosity = 1.81e-5,       # Pa·s at STP
            thermal_conductivity = 0.025,  # W/m·K at STP
            specific_heat = 1007.0,    # J/kg·K
            thermal_expansion = 3.67e-3,   # 1/K at STP
            compressibility = 1e-6     # 1/Pa (nearly incompressible)
        )
        
    elseif model == :temperature_dependent
        # Sutherland's law for air viscosity
        viscosity_model = SutherlandModel(
            reference_viscosity = 1.716e-5,      # Pa·s at 273K
            reference_temperature = 273.0,       # K
            sutherland_constant = 110.4          # K for air
        )
        
        # Thermal conductivity correlation for air
        thermal_conductivity_model = PowerLawModel(
            coefficient = 2.334e-3,    # W/m·K·K^(-0.7)
            exponent = 0.7             # Temperature exponent
        )
        
        # Ideal gas law for density
        density_model = IdealGasModel(
            molar_mass = 0.02897,      # kg/mol (dry air)
            gas_constant = 8.314       # J/mol·K
        )
        
        # Specific heat correlation for air
        specific_heat_model = PolynomialModel(
            coefficients = [1047.63, -0.372589, 9.45304e-4, -6.02409e-7]  # cₚ(T) polynomial
        )
        
        return TemperatureDependentProperties(
            viscosity_model = viscosity_model,
            thermal_conductivity_model = thermal_conductivity_model,
            density_model = density_model,
            specific_heat_model = specific_heat_model,
            reference_conditions = (
                T₀ = 273.15,    # K
                P₀ = 101325.0,  # Pa
                ρ₀ = 1.292,     # kg/m³
                μ₀ = 1.716e-5,  # Pa·s
                k₀ = 0.0241,    # W/m·K
                cₚ₀ = 1005.0    # J/kg·K
            ),
            temperature_range = (200.0, 2000.0),
            pressure_range = (1e3, 1e7)
        )
        
    elseif model == :ideal_gas
        # Simple ideal gas properties
        return TemperatureDependentProperties(
            viscosity_model = SutherlandModel(1.716e-5, 273.0, 110.4),
            thermal_conductivity_model = PowerLawModel(2.334e-3, 0.7),
            density_model = IdealGasModel(),
            specific_heat_model = PowerLawModel(1007.0, 0.0),  # Constant cₚ
            reference_conditions = (T₀=273.15, P₀=101325.0, ρ₀=1.292, μ₀=1.716e-5),
            temperature_range = (200.0, 1500.0),
            pressure_range = (1e3, 1e6)
        )
    else
        throw(ArgumentError("Unknown air property model: $model"))
    end
end

"""
    create_water_properties(; model=:temperature_dependent) -> FluidProperties

Create water property database with temperature-dependent correlations.

# Arguments
- `model::Symbol`: Property model type (:constant, :temperature_dependent)

# Returns
- `FluidProperties`: Water properties for liquid water
"""
function create_water_properties(; model::Symbol=:temperature_dependent)
    if model == :constant
        return ConstantProperties(
            density = 998.2,           # kg/m³ at 20°C
            viscosity = 1.002e-3,      # Pa·s at 20°C
            thermal_conductivity = 0.598,  # W/m·K at 20°C
            specific_heat = 4182.0,    # J/kg·K
            thermal_expansion = 2.14e-4,   # 1/K at 20°C
            compressibility = 4.6e-10  # 1/Pa (liquid water)
        )
        
    elseif model == :temperature_dependent
        # Temperature-dependent water properties
        
        # Viscosity: exponential model μ(T) = A*exp(B/T)
        viscosity_model = CustomFunctionModel(
            property_function = (T, P=1e5) -> 2.414e-5 * 10^(247.8/(T-140)),
            parameter_dict = Dict("type" => "water_viscosity", "range" => (273.15, 373.15))
        )
        
        # Thermal conductivity: polynomial model
        thermal_conductivity_model = PolynomialModel(
            coefficients = [-0.432, 0.0057, -6.4e-6, 1.7e-9]  # k(T) in W/m·K
        )
        
        # Density: polynomial model for liquid water
        density_model = PolynomialModel(
            coefficients = [1002.56, -0.1262, -3.79e-3, 2.05e-6]  # ρ(T) in kg/m³
        )
        
        # Specific heat: polynomial model
        specific_heat_model = PolynomialModel(
            coefficients = [4217.4, -3.720616, 0.0138, -2.24e-5]  # cₚ(T) in J/kg·K
        )
        
        return TemperatureDependentProperties(
            viscosity_model = viscosity_model,
            thermal_conductivity_model = thermal_conductivity_model,
            density_model = density_model,
            specific_heat_model = specific_heat_model,
            reference_conditions = (
                T₀ = 293.15,    # K (20°C)
                P₀ = 101325.0,  # Pa
                ρ₀ = 998.2,     # kg/m³
                μ₀ = 1.002e-3,  # Pa·s
                k₀ = 0.598,     # W/m·K
                cₚ₀ = 4182.0    # J/kg·K
            ),
            temperature_range = (273.15, 373.15),  # Liquid water range
            pressure_range = (1e3, 1e6)
        )
    else
        throw(ArgumentError("Unknown water property model: $model"))
    end
end

"""
    create_steam_properties() -> FluidProperties

Create steam (water vapor) property database.
"""
function create_steam_properties()
    # Steam as ideal gas with temperature-dependent properties
    
    # Viscosity from kinetic theory
    viscosity_model = SutherlandModel(
        reference_viscosity = 1.34e-5,    # Pa·s at 373K
        reference_temperature = 373.0,    # K
        sutherland_constant = 1064.0      # K for steam
    )
    
    # Thermal conductivity
    thermal_conductivity_model = PowerLawModel(
        coefficient = 1.614e-2,   # Coefficient for steam
        exponent = 1.0            # Linear with temperature
    )
    
    # Ideal gas density
    density_model = IdealGasModel(
        molar_mass = 0.018015,    # kg/mol (H₂O)
        gas_constant = 8.314      # J/mol·K
    )
    
    # Specific heat (temperature-dependent)
    specific_heat_model = PolynomialModel(
        coefficients = [1863.0, 0.441, -1.5e-4]  # Steam cₚ(T)
    )
    
    return TemperatureDependentProperties(
        viscosity_model = viscosity_model,
        thermal_conductivity_model = thermal_conductivity_model,
        density_model = density_model,
        specific_heat_model = specific_heat_model,
        reference_conditions = (
            T₀ = 373.15,     # K
            P₀ = 101325.0,   # Pa
            ρ₀ = 0.598,      # kg/m³
            μ₀ = 1.34e-5,    # Pa·s
            k₀ = 0.025,      # W/m·K
            cₚ₀ = 2010.0     # J/kg·K
        ),
        temperature_range = (373.15, 1273.15),  # Steam range
        pressure_range = (1e3, 1e7)
    )
end

# Property evaluation functions

"""
    compute_viscosity(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute dynamic viscosity at given temperature and pressure.

# Arguments
- `T::Float64`: Temperature (K)
- `P::Float64`: Pressure (Pa)
- `properties::FluidProperties`: Fluid property database

# Returns
- `Float64`: Dynamic viscosity μ (Pa·s)
"""
function compute_viscosity(T::Float64, P::Float64, properties::FluidProperties)
    if properties isa ConstantProperties
        return properties.viscosity
        
    elseif properties isa TemperatureDependentProperties
        return evaluate_property_model(T, P, properties.viscosity_model)
        
    elseif properties isa NonNewtonianProperties
        # Base viscosity from underlying model
        μ₀ = compute_viscosity(T, P, properties.base_properties)
        
        # Apply non-Newtonian correction (requires shear rate)
        # For now, return base viscosity - full non-Newtonian would need γ̇
        return μ₀
    else
        throw(ArgumentError("Unknown fluid properties type: $(typeof(properties))"))
    end
end

"""
    compute_thermal_conductivity(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute thermal conductivity at given temperature and pressure.
"""
function compute_thermal_conductivity(T::Float64, P::Float64, properties::FluidProperties)
    if properties isa ConstantProperties
        return properties.thermal_conductivity
        
    elseif properties isa TemperatureDependentProperties
        return evaluate_property_model(T, P, properties.thermal_conductivity_model)
        
    elseif properties isa NonNewtonianProperties
        return compute_thermal_conductivity(T, P, properties.base_properties)
        
    else
        throw(ArgumentError("Unknown fluid properties type: $(typeof(properties))"))
    end
end

"""
    compute_density(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute density at given temperature and pressure.
"""
function compute_density(T::Float64, P::Float64, properties::FluidProperties)
    if properties isa ConstantProperties
        return properties.density
        
    elseif properties isa TemperatureDependentProperties
        return evaluate_property_model(T, P, properties.density_model)
        
    elseif properties isa NonNewtonianProperties
        return compute_density(T, P, properties.base_properties)
        
    else
        throw(ArgumentError("Unknown fluid properties type: $(typeof(properties))"))
    end
end

"""
    compute_specific_heat(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute specific heat capacity at given temperature and pressure.
"""
function compute_specific_heat(T::Float64, P::Float64, properties::FluidProperties)
    if properties isa ConstantProperties
        return properties.specific_heat
        
    elseif properties isa TemperatureDependentProperties
        return evaluate_property_model(T, P, properties.specific_heat_model)
        
    elseif properties isa NonNewtonianProperties
        return compute_specific_heat(T, P, properties.base_properties)
        
    else
        throw(ArgumentError("Unknown fluid properties type: $(typeof(properties))"))
    end
end

"""
    compute_thermal_diffusivity(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute thermal diffusivity α = k/(ρ·cₚ).
"""
function compute_thermal_diffusivity(T::Float64, P::Float64, properties::FluidProperties)
    k = compute_thermal_conductivity(T, P, properties)
    ρ = compute_density(T, P, properties)
    cₚ = compute_specific_heat(T, P, properties)
    
    return k / (ρ * cₚ)
end

"""
    compute_kinematic_viscosity(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute kinematic viscosity ν = μ/ρ.
"""
function compute_kinematic_viscosity(T::Float64, P::Float64, properties::FluidProperties)
    μ = compute_viscosity(T, P, properties)
    ρ = compute_density(T, P, properties)
    
    return μ / ρ
end

"""
    compute_prandtl_number(T::Float64, P::Float64, properties::FluidProperties) -> Float64

Compute Prandtl number Pr = ν/α = (μ·cₚ)/k.
"""
function compute_prandtl_number(T::Float64, P::Float64, properties::FluidProperties)
    μ = compute_viscosity(T, P, properties)
    k = compute_thermal_conductivity(T, P, properties)
    cₚ = compute_specific_heat(T, P, properties)
    
    return (μ * cₚ) / k
end

# Property model evaluation

"""
    evaluate_property_model(T::Float64, P::Float64, model::PropertyModel) -> Float64

Evaluate a property correlation model at given conditions.
"""
function evaluate_property_model(T::Float64, P::Float64, model::PropertyModel)
    if model isa SutherlandModel
        μ₀ = model.reference_viscosity
        T₀ = model.reference_temperature
        S = model.sutherland_constant
        
        return μ₀ * (T/T₀)^1.5 * (T₀ + S)/(T + S)
        
    elseif model isa PowerLawModel
        A = model.coefficient
        n = model.exponent
        
        return A * T^n
        
    elseif model isa PolynomialModel
        coeffs = model.coefficients
        result = 0.0
        
        for (i, coeff) in enumerate(coeffs)
            result += coeff * T^(i-1)
        end
        
        return result
        
    elseif model isa PiecewiseLinearModel
        T_points = model.temperature_points
        prop_values = model.property_values
        
        # Linear interpolation
        if T <= T_points[1]
            return prop_values[1]
        elseif T >= T_points[end]
            return prop_values[end]
        else
            # Find interpolation interval
            for i in 1:(length(T_points)-1)
                if T_points[i] <= T <= T_points[i+1]
                    # Linear interpolation
                    t = (T - T_points[i]) / (T_points[i+1] - T_points[i])
                    return prop_values[i] + t * (prop_values[i+1] - prop_values[i])
                end
            end
        end
        
    elseif model isa IdealGasModel
        M = model.molar_mass
        R = model.gas_constant
        
        return P * M / (R * T)
        
    elseif model isa CustomFunctionModel
        return model.property_function(T, P)
        
    else
        throw(ArgumentError("Unknown property model: $(typeof(model))"))
    end
end

# Non-Newtonian rheology

"""
    compute_apparent_viscosity(shear_rate::Float64, T::Float64, properties::NonNewtonianProperties) -> Float64

Compute apparent viscosity for non-Newtonian fluids.

# Arguments
- `shear_rate::Float64`: Magnitude of shear rate tensor γ̇ (1/s)
- `T::Float64`: Temperature (K)
- `properties::NonNewtonianProperties`: Non-Newtonian properties

# Returns
- `Float64`: Apparent viscosity μₐ (Pa·s)
"""
function compute_apparent_viscosity(shear_rate::Float64, T::Float64, properties::NonNewtonianProperties)
    model = properties.rheology_model
    params = properties.rheology_parameters
    
    # Base viscosity (Newtonian component)
    μ₀ = compute_viscosity(T, 101325.0, properties.base_properties)
    
    if model == :power_law
        # Power law: μₐ = K * γ̇^(n-1)
        K = get(params, :consistency_index, μ₀)  # Consistency index
        n = get(params, :flow_behavior_index, 1.0)  # Flow behavior index
        
        if abs(shear_rate) < 1e-12
            return K * 1e12^(n-1)  # Large viscosity for zero shear rate
        else
            return K * abs(shear_rate)^(n-1)
        end
        
    elseif model == :bingham
        # Bingham plastic: μₐ = μₚ + τᵧ/γ̇
        μₚ = get(params, :plastic_viscosity, μ₀)  # Plastic viscosity
        τᵧ = get(params, :yield_stress, 0.0)      # Yield stress
        
        if abs(shear_rate) < 1e-12
            return 1e12  # Infinite viscosity below yield
        else
            return μₚ + τᵧ / abs(shear_rate)
        end
        
    elseif model == :herschel_bulkley
        # Herschel-Bulkley: μₐ = (τᵧ/γ̇) + K*γ̇^(n-1)
        τᵧ = get(params, :yield_stress, 0.0)
        K = get(params, :consistency_index, μ₀)
        n = get(params, :flow_behavior_index, 1.0)
        
        if abs(shear_rate) < 1e-12
            return 1e12
        else
            return τᵧ / abs(shear_rate) + K * abs(shear_rate)^(n-1)
        end
        
    else
        @warn "Unknown rheology model: $model, using Newtonian"
        return μ₀
    end
end

# Mixture properties

"""
    compute_mixture_properties(mass_fractions::Vector{Float64}, component_properties::Vector{FluidProperties},
                              T::Float64, P::Float64) -> NamedTuple

Compute mixture properties from component properties and mass fractions.

# Arguments
- `mass_fractions::Vector{Float64}`: Mass fractions of components (must sum to 1)
- `component_properties::Vector{FluidProperties}`: Properties of each component
- `T::Float64`: Temperature (K)
- `P::Float64`: Pressure (Pa)

# Returns
- `NamedTuple`: Mixture properties (density, viscosity, thermal_conductivity, specific_heat)
"""
function compute_mixture_properties(mass_fractions::Vector{Float64}, 
                                   component_properties::Vector{FluidProperties},
                                   T::Float64, P::Float64)
    if abs(sum(mass_fractions) - 1.0) > 1e-6
        throw(ArgumentError("Mass fractions must sum to 1.0"))
    end
    
    if length(mass_fractions) != length(component_properties)
        throw(ArgumentError("Number of mass fractions must match number of components"))
    end
    
    n_components = length(mass_fractions)
    
    # Compute component properties
    ρᵢ = [compute_density(T, P, props) for props in component_properties]
    μᵢ = [compute_viscosity(T, P, props) for props in component_properties]
    kᵢ = [compute_thermal_conductivity(T, P, props) for props in component_properties]
    cₚᵢ = [compute_specific_heat(T, P, props) for props in component_properties]
    
    # Mixture density (mass-weighted harmonic mean)
    ρ_mix = 1.0 / sum(mass_fractions[i] / ρᵢ[i] for i in 1:n_components)
    
    # Mixture viscosity (mass-weighted)
    μ_mix = sum(mass_fractions[i] * μᵢ[i] for i in 1:n_components)
    
    # Mixture thermal conductivity (mass-weighted)
    k_mix = sum(mass_fractions[i] * kᵢ[i] for i in 1:n_components)
    
    # Mixture specific heat (mass-weighted)
    cₚ_mix = sum(mass_fractions[i] * cₚᵢ[i] for i in 1:n_components)
    
    return (
        density = ρ_mix,
        viscosity = μ_mix,
        thermal_conductivity = k_mix,
        specific_heat = cₚ_mix,
        thermal_diffusivity = k_mix / (ρ_mix * cₚ_mix),
        kinematic_viscosity = μ_mix / ρ_mix,
        prandtl_number = (μ_mix * cₚ_mix) / k_mix
    )
end

# Property validation and bounds checking

"""
    validate_properties(T::Float64, P::Float64, properties::FluidProperties) -> Bool

Validate that computed properties are physically reasonable.

# Arguments
- `T::Float64`: Temperature (K)
- `P::Float64`: Pressure (Pa)
- `properties::FluidProperties`: Fluid properties to validate

# Returns
- `Bool`: true if properties are valid, false otherwise
"""
function validate_properties(T::Float64, P::Float64, properties::FluidProperties)
    # Check temperature and pressure ranges
    if properties isa TemperatureDependentProperties
        T_min, T_max = properties.temperature_range
        P_min, P_max = properties.pressure_range
        
        if T < T_min || T > T_max
            @warn "Temperature $T K outside valid range [$T_min, $T_max] K"
            return false
        end
        
        if P < P_min || P > P_max
            @warn "Pressure $P Pa outside valid range [$P_min, $P_max] Pa"
            return false
        end
    end
    
    # Compute properties and check for physical reasonableness
    try
        ρ = compute_density(T, P, properties)
        μ = compute_viscosity(T, P, properties)
        k = compute_thermal_conductivity(T, P, properties)
        cₚ = compute_specific_heat(T, P, properties)
        
        # Physical bounds checking
        if ρ <= 0 || ρ > 20000  # Density bounds (kg/m³)
            @warn "Unphysical density: $ρ kg/m³"
            return false
        end
        
        if μ <= 0 || μ > 1000  # Viscosity bounds (Pa·s)
            @warn "Unphysical viscosity: $μ Pa·s"
            return false
        end
        
        if k <= 0 || k > 1000  # Thermal conductivity bounds (W/m·K)
            @warn "Unphysical thermal conductivity: $k W/m·K"
            return false
        end
        
        if cₚ <= 0 || cₚ > 50000  # Specific heat bounds (J/kg·K)
            @warn "Unphysical specific heat: $cₚ J/kg·K"
            return false
        end
        
        return true
        
    catch e
        @warn "Error computing properties: $e"
        return false
    end
end

# Utility functions

"""
    create_custom_fluid(; density_model, viscosity_model, thermal_conductivity_model, 
                        specific_heat_model, reference_conditions, 
                        temperature_range, pressure_range) -> TemperatureDependentProperties

Create custom fluid properties with user-defined correlations.
"""
function create_custom_fluid(; density_model::PropertyModel,
                            viscosity_model::PropertyModel,
                            thermal_conductivity_model::PropertyModel,
                            specific_heat_model::PropertyModel,
                            reference_conditions::NamedTuple,
                            temperature_range::Tuple{Float64,Float64} = (200.0, 2000.0),
                            pressure_range::Tuple{Float64,Float64} = (1e3, 1e7))
    
    return TemperatureDependentProperties(
        density_model = density_model,
        viscosity_model = viscosity_model,
        thermal_conductivity_model = thermal_conductivity_model,
        specific_heat_model = specific_heat_model,
        reference_conditions = reference_conditions,
        temperature_range = temperature_range,
        pressure_range = pressure_range
    )
end

"""
    print_property_summary(T::Float64, P::Float64, properties::FluidProperties)

Print a summary of computed fluid properties at given conditions.
"""
function print_property_summary(T::Float64, P::Float64, properties::FluidProperties)
    println("Fluid Properties Summary")
    println("========================")
    println("Temperature: $(T) K ($(T-273.15)°C)")
    println("Pressure: $(P) Pa ($(P/1e5) bar)")
    println()
    
    if validate_properties(T, P, properties)
        ρ = compute_density(T, P, properties)
        μ = compute_viscosity(T, P, properties)
        ν = compute_kinematic_viscosity(T, P, properties)
        k = compute_thermal_conductivity(T, P, properties)
        cₚ = compute_specific_heat(T, P, properties)
        α = compute_thermal_diffusivity(T, P, properties)
        Pr = compute_prandtl_number(T, P, properties)
        
        println("Density (ρ):                 $(ρ) kg/m³")
        println("Dynamic viscosity (μ):       $(μ) Pa·s")
        println("Kinematic viscosity (ν):     $(ν) m²/s")
        println("Thermal conductivity (k):    $(k) W/m·K")
        println("Specific heat (cₚ):          $(cₚ) J/kg·K")
        println("Thermal diffusivity (α):     $(α) m²/s")
        println("Prandtl number (Pr):         $(Pr)")
    else
        println("ERROR: Invalid properties at these conditions")
    end
end